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We  developed  a  spiCE-compatible  compact  model  of  TiCE-TiCE  ,  memristors  based  on  classic  ion 
transportation  theory.  Our  model  is  shown  to  simulate  important  dynamic  memristive  properties 
like  real-time  memristance  switching,  which  are  critical  in  memristor-based  analog  circuit  designs. 
The  model,  as  well  as  its  analytical  approximation,  is  validated  with  the  experimentally  obtained 
data  from  real  devices.  Minor  deviations  of  our  model  from  the  measured  data  are  also  analyzed 
and  discussed.  ©2013  AIP  Publishing  LLC  [http://dx.doi.org/10.1063/L4802206] 


The  operational  characteristics  of  memristive  devices 
are  typically  represented  by  an  I-V  curve:  the  memristance 
of  the  device  changes  with  the  magnitude  and  pulse  duration 
of  the  externally  applied  excitations.1  As  one  of  the  promis¬ 
ing  technologies,  Ti02-Ti02_A  memristor  has  recently 
received  significant  attention  and  been  widely  studied  in  the 
solid  state  device  society. 2-10  Some  physical  models  of 
memristive  devices,  such  as  HfO,  and  NiO,,  proposed  a  gen¬ 
eral  description  on  the  switching  on/off  procedure  based  on 
the  ion/vacancy  motion  driven  by  the  electric  field  (or  poten¬ 
tial  gradient),  which  is  known  as  the  filament  formation  or 
dissolution. 11-18  The  behavior  of  the  TiCE-TiCE-*  device  can 
also  be  understood  by  the  similar  ionic  influence.  ’  the 
filament  corresponds  to  the  high  conductive  region  while  the 
remaining  regions  can  be  treated  as  low  conductive  or  insu¬ 
lating.  Some  other  models  present  the  understandings  on  in¬ 
ternal  state  changes  under  physical  factors  like  electric  field 
and  temperature.1 112  24  Some  previous  work  has  proposed  to 
model  the  memristive  switching  behavior  in  oxide  devices 
based  on  the  dynamical  interactions  between  electrons,  oxy¬ 
gen  vacancies,  and  oxygen  anions.22  However,  due  to  the 
high  computation  cost  of  complex  physical  quantities,  i.e., 
transition  Hamilton  matrix  and  energy  density  of  states,  the 
model  cannot  be  directly  applicable  to  fast  circuit  simulation. 
By  combining  Monte-Carlo  method,  electron  hopping  con¬ 
duction  process  can  be  simulated  through  the  similar  tech¬ 
nique.26  Kinetic  Monte-Carlo  method  is  also  implemented  to 
model  the  filament  formation.27  However,  these  models  do 
not  explicitly  give  the  static  and  dynamic  electrical  switch¬ 
ing  properties  of  the  memristor  device.  In  some  TiOx  mem¬ 
ristor  modeling  works,  the  motion  of  vacancies  is  simulated 
through  classical  molecular  dynamics.28  29  The  dynamics  of 
tunnel  barrier  width  can  be  also  obtained  by  modeling  the 
device  as  a  resistor  plus  Simmons  tunnel  barrier.30  In  this 
work,  we  model  the  behavior  of  memristor  device  from  clas¬ 
sical  macroscopic  viewpoint,  which  is  represented  by  the 
electron  density  distribution  inside  the  device.  A  compact 
model  of  TiCE-TiCE-*  memristor  based  on  classic  ion  trans¬ 
portation  theory  is  then  proposed,  targeting  the  application  in 


spice  circuit  simulations.  Compared  to  the  aforementioned 
TiOv  memristor  models,  our  model  offers  a  concrete  calcula¬ 
tion  of  the  filament  dynamics  and  derives  accurate  time- 
varying  behaviors  of  device  resistance  changing.  Our  model, 
as  well  as  its  analytical  approximation,  well  matches  the 
measurements  of  a  real  TiCE-TiCE-*  device  on  static  I-V 
curve  and  dynamic  pulse  programming. 

We  divide  a  Ti02-Ti02-x  memristor  into  three  regions, 
namely,  the  conductive  region,  the  transition  region,  and  the 
insulating  region,  as  shown  in  Fig.  1.  Here,  the  conductive 
region  corresponds  to  the  region  of  filaments  while  the  tran¬ 
sition  region  connects  the  filaments  and  the  remaining  insu¬ 
lating  region.  The  formation  or  dissolution  of  the  filament  is 
simulated  as  the  ion  surface  motion  within  the  transition 
region.  We  use  w,  1,  and  D  to  denote  the  lengths  of  the  con¬ 
ductive  region,  the  transition  region,  and  the  entire  device, 
respectively.  The  real-time  current  density  at  position  x  in 
the  memristor  at  time  t  can  be  expressed  as3 1 

,  „  ,  ,  ,  ,  dnJx.t) 

Js(x,  t )  =  ns(x,  t)qpEs(x,  t )  -  qDq  — — — .  (1) 

Here,  the  subscript  s  =  ( c ,  t,  i )  denotes  the  parameter  of  con¬ 
ductive  region,  transition  region,  and  insulating  region, 


V(t) 
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FIG.  1.  Region  partitioning  of  TiCE-TiCE-A- memristor. 
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respectively.  Es  and  ns  represent  the  electric  field  and  the 
electron  density,  respectively,  fi  is  the  mobility  coefficient,  q 
is  the  elementary  charge.  Dq  is  the  diffusion  constant  for 
electrons.  Equation  (1)  shows  that  the  current  is  generated 
mainly  from  the  electron  drifting  in  the  electric  field  and  the 
electron  density  gradient.  Dq  is  typically  small  for  electron 
diffusion  current,  and  the  second  term  in  Eq.  (1)  is  ignored 
in  our  model.  If  we  ignore  the  variation  of  the  cross  section 
area  and  assume  the  current  density  is  uniform  in  the  mem- 
ristor  (i.e.,  Jc  =  J,  =  /,),  the  relationships  between  the  elec¬ 
tric  fields  and  electron  densities  of  the  three  regions  can  be 
summarized  as 

ncqfiEc(t )  =  nt(x,t)qfiEt(x,t)  =  ni{t)qfiEj{t).  (2) 

Here,  we  assume  that  the  Es  and  ns  are  uniform  in  conductive 
and  insulating  regions  and  are  determined  only  by  the 
real-time  applied  voltage  (except  that  nc  is  a  constant).  The 
equation  for  the  hopping  rate  of  ions  in  a  infinitely  large  lat¬ 
tice  can  be  expressed  as32 

R  =  fsmh{QEa/2kT)e~u/kT,  (3) 

where  R  is  the  ion  hopping  rate, /is  the  escape  attempt  fre¬ 
quency,  a  is  the  lattice  periodicity,  E  is  the  electric  field,  Q  is 
the  charge  of  ion,  U  is  the  activation  energy,  k  is  the 
Boltzmann  constant,  and  T  is  the  local  temperature.  By  con¬ 
sidering  the  fact  that  the  number  of  generated  electrons  is 
proportional  to  that  of  the  hopping  ions  inside  the  thin  film, 
we  use  the  product  of  a  coefficient  y,  and  a  factor  (nc  —  «,(r)) 
to  approximate  the  electron  density  generation  function  under 
the  influence  of  electric  field  in  the  insulating  region.  ( nc 
—rii(t))  also  implies  that  the  maximum  electron  density 
should  not  exceed  nc.  Hence,  the  evolution  of  the  electron 
density  in  insulating  region  under  a  time-varying  electric  field 
Ej  can  be  calculated  by 

^^  =  7/K-«/(0)£/(0,  (4) 

where  yi  is  the  electron  density  generating  coefficient  in  the 
insulating  region.  At  a  given  time  t  and  position  x,  the  elec¬ 
tric  field  Et  can  be  viewed  as  a  linear  function  bounded  by  Er 
and  Ej  as 


Et(x,  t ) 


(£,-(Q-£c(Q) 

A 


(x-w)  +Ec(t). 


(5) 


We  note  that  the  voltage  V(t)  applied  on  the  two  ends  of 
the  memristor  equals  the  integral  of  the  electric  field  along 
the  device,  or 


W+A 


V(t)  =  Ec(t)w  - 


Et(x,  t)dx  +  Ej(t)(D  —  A  —  w).  (6) 


Combining  Eqs.  (5)  and  (6),  we  have 

V{i) 


Ec(t)  = 


1  nr 

W  +  -  A  -r— 

2  m{t) 


D  —  w  —  -A 


(7) 


Finally,  the  transition  region  length  A  reduces  when  the 
applied  voltage  V(t )  increases,  which  is  approximated  by 

A  =  Aoe~  |K(,)I,  (8) 


in  our  model.  Aq  is  the  transition  region  length  at  V(t)  =  0. 

The  probability  of  ion  hopping  changes  with  the  external 
electric  field.12  However,  the  probability  is  significantly  dif¬ 
ferent  in  the  conductive  and  insulating  regions,  which  lead¬ 
ing  to  the  conductivity  difference  in  these  two  regions. 
Ideally,  every  successful  ion  hopping  produces  two  free  elec¬ 
trons  and  can  be  viewed  as  the  electron  density  increase  or 
redistribution.  For  simplicity,  we  assume  the  conductivity  in 
the  conductive  region  (filament)  is  constant.  The  motion  of 
ions  causes  the  change  of  the  electron  densities  in  transition 
and  insulating  regions,  resulting  in  the  formation,  growth, 
and  dissolution  of  the  filaments.  Since  the  TiCE-TiCE_x 
device  is  stable  when  no  voltage  bias  is  applied,  we  can 
model  the  ion  motion  by  considering  only  the  external  elec¬ 
tric  field.12  The  internal  temperature  of  the  device,  which 
leads  to  the  fluctuations  of  cycle-to-cycle  switching,  is  not 
included  in  our  model.  Without  loss  of  generality,  we  assume 
the  endpoint  of  the  filament  starts  moving  from  the  left  end 
of  the  device  (x=  0)  at  time  0.  nt{x,  t)  is  the  electron  density 
at  the  position  x  in  the  transition  region  at  time  t.  Compared 
to  the  electron  density  nc  at  the  boundary  between  the  con¬ 
ductive  region  and  the  transition  region,  the  change  of  the 
electron  density  at  jc  is  nc  —  n,(x.  t),  which  is  mainly  due  to 
the  oxygen  ion  vacancy  redistribution.33 

When  external  voltage  is  applied,  the  electron  densities 
are  changing  with  time  in  both  transition  and  insulating 
regions  (filament  growth).  The  electron  density  increment 
in  an  infinitesimal  time  interval  dt  equals  the  difference 
between  the  electron  densities  at  positions  x  -  dx  and  x,  or 


,  .  .  .  dn,( x,  t)  , 

nt(x  —  dx,  t)  —  nt(x ,  t )  =  — — - {—dx) 


dx 
dn,(x ,  t) 
dx 


v{x,  t)dt 


=  y,{nc  -  n,{x,  t))E,{x ,  t)dt.  (9) 


Here,  y,  is  the  electron  generating  coefficient  in  the  transition 
region.  Note  that  the  maximum  reachable  electron  density  in 
the  memristor  device  is  nc,  or  the  region  is  fully  conductive. 

Based  on  Eq.  (9),  the  growth  velocity  at  position  x  can 
be  calculated  by 


y,{nc  -  n,(x,  t))E,(x ,  t) 
dnt{x,  t) 
dx 


(10) 


where  d"'T’t’  can  be  derived  from  Eqs.  (2)  and  (5). 

The  motion  of  the  transition  region  can  be  described  by 
the  average  growth  velocity,  which  is  defined  as 

w-\-A 

v(x,  t)dx 

v{t)  = 


A 


(11) 
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This  result  shows  that  the  filament  growth  is  triggered  by 
mainly  the  external  electric  field  and  can  be  viewed  as  the 
combined  imaginary  surfaces  within  the  transition  region.  In 
the  following  sections,  we  simplify  the  expressions  of  nc, 
nt(x,  t ),  and  n,(?)  as  nc,  nt,  and  which  are  still  the  functions 
of  x  and/or  t. 

Substituting  Eqs.  (2)  and  (5)  into  Eq.  (10),  we  have 


v(x,  t) 


ncEc(t)  ■ 


7,  knciii 

nj  nc  -  tii ' 


(12) 


Substituting  Eqs.  (2),  (5),  (7),  and  ( 12)  into  Eq.  (11),  the  fila¬ 
ment  growth  velocity  can  be  approximated  by 


dw(t) 

dt 


v(t) 


=  7tW  ■ 


«/ 


V(t) 


U+- 

2  fi; 


(13) 


Here,  /I  equals  I  (|  -  l)3  +  \  (*  -  l)2  +  \  (|  -  1). 

The  resistance/memristance  of  the  memristor  can  be  cal¬ 
culated  by 


*(0  = 


JCA 


V{t) 

ncqfiEc(t)A 


(14) 


Here,  A  is  the  cross  section  area  of  the  memristor  filament. 
Equations  (13)  and  (14)  describe  the  dynamic  changes  of  mem¬ 
ristor  device  structure  and  electrical  property,  respectively. 

Assuming  nc  and  n,  are  constant  during  the  switching  of 
the  memristor  (for  example,  the  applied  voltage  magnitude  is 
low  or  the  programming  pulse  width  is  short),  we  can  obtain 
the  general  form  of  w  under  a  constant  applied  voltage  V  by 
integrating  both  sides  of  Eq.  (13)  as 


W(t)  : 


1  — 


2 


w(fo)  +  Co 


2yt[l/,V At  —  Co  }. 


(15) 


Here,  At  =  t  —  to  and  Co  =  \  2  +  J  (D  —  1 2).  Equation  (15) 
provides  an  analytical  approximation  on  the  growth  of  the  fil¬ 
aments  in  the  memristor.  The  internal  state  variable  w  is  the 
expression  of  electron  densities  and  the  previous  value  of  w0. 
This  result  bridges  two  state  variables  at  two  different  times 
under  one  fixed  voltage,  w  represents  the  filament  length,  or 
the  state  of  the  device  and  is  able  to  calculate  the  state  vari¬ 
able  value  or  the  total  resistance.  As  we  shall  see  later,  this 
analytical  result  is  close  to  the  experimental  result  of  the  total 
device  resistance.  Based  on  that,  we  can  quickly  estimate  the 
memristor  switching  time  in  the  following  two  scenarios: 

If  the  device  is  switched  on,  or  the  filament  is  formed 
from  the  left-end  of  the  memristor  to  the  right-end,  at  ?o  =  0 
(or  w(?o)  =  0),  Eq.  (15)  can  be  simplified  as 


1  -  — 


The  memristor  switching  time  can  be  estimated  by  assuming 
Cl  -  2 ytfVkt  =  0,  or 


ton  — 


c2 


2y,m 


(17) 


The  corresponding  portion  of  transition  region  at  r0  is  w(ton) 
k.D  —  ^2  since  n,  <C  nc.  Considering  normally  /.  <C  IK  it 
indicates  that  the  filament  is  formed  completely  at  the  right 
end  of  the  memristor  at  t  =  t„„ ,  or  the  device  is  totally 
switched  on. 

Similarly,  if  the  filament  is  in  the  dissolution  process 
starts  from  the  right-end  of  the  memristor  at  to  =  0  (or  w(0) 
=  D),  Eq.  (15)  can  be  simplified  as 


w{t)  = 


1 


2y, fi/.Vt  -  C 


(18) 


The  time  required  to  switch  off  the  memristor  (t0g)  can  be 
derived  by  assuming  h '(toff)  =  0,  or 


toff  = 


(  nc\ 

(  nc\ 

1  --  D  +  2C0 

( 1  — -)d 

LV  «.7 

2  ytp2V 


(19) 


Table  I  summarized  the  three  types  of  the  device  param¬ 
eters  used  in  our  model,  including  the  geometric  parameters, 
the  electrical  parameters,  and  the  structural  parameters.  The 
electron  generating  coefficients  y,  and  y,  are  derived  from 
the  measured  data  and  assumed  constant  for  the  different 
working  ranges,  as  shown  in  Fig.  2.  We  compared  our  model 
with  the  experimentally  obtained  characterized  static  I-V 
curve  and  dynamic  pulse  programming  curve  of  a  Ti02- 
Ti02_A  memristor  device. 

Fig.  2  shows  the  measured  I-V  curve  from  a  real  mem¬ 
ristor  device  as  well  as  the  experimental  data  from  the  nu¬ 
merical  simulation  and  analytical  approximation.  During  the 
measurement,  a  sequence  of  voltage  pulses  is  applied  on  the 
memristor  device.  The  magnitude  of  the  voltage  pulse  grows 
exponentially,  and  varies  from  positive  to  negative  following 
sine  function.  The  numerical  simulation  well  fit  the  meas¬ 
ured  data  in  all  four  working  ranges.  The  analytical  approxi¬ 
mation  shows  slight  discrepancy  from  the  measured  data 
when  the  resistance  is  high.  It  is  because  the  simulated 


TABLE  I.  Modeling  parameters. 


Parameters 

Value 

Parameters 

Value 

Geometric 

D 

35  nm 

A 

0.25  /(m2 

Electrical 

e 

1.602  x  10“19C 

nc 

8.75  x  1022  m~3 

Structural 

w0 

0.15  D 

0.05  D 

Working  range 

Derived  parameters 

1  ->  2 

y, 

2.3  x  10“6 

Vi 

1  x  10~6 

2  -t-  1 

8  x  10“6 

1  x  10'8 

1  3 

1  x  10-10 

1  x  10~9 

3  -►  1 

7  x  10^7 

2  x  10~7 

w(t)  = 


Cq  —  2ytflVXt  —  C0 


(16) 
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FIG.  2.  Model  validation  with  static  I-V  curve,  including  numerical  simula¬ 
tion  and  analytical  approximation.  The  magnitude  of  the  applied  voltage 
grows  exponentially  and  follows  sine  function.  The  memristor  is  initially  set 
to  high  resistance  state. 

filament  growth  velocity  is  lower  than  the  actual  value  when 
the  variation  of  n,  over  time  is  ignored. 

To  prove  the  capability  of  our  model  on  simulating  the 
dynamic  switching  property  of  the  memristor,  we  plot  the  re¬ 
sistance  changes  following  the  programming  pulses  in 
Fig.  3.  The  resistance  of  the  memristor  first  decreases  when 
the  positive  pulses  are  applied,  and  then  raises  when  the 
polarization  of  the  pulses  changes  to  negative.  Our  numerical 
simulation  matches  the  measured  data  very  well  over  most 
of  the  plotted  points.  Small  discrepancies  show  at  the  high 
resistance  state.  One  reason  for  the  deviations  could  be  the 
impact  of  thermal  fluctuations,  which  become  prominent 
under  a  relatively  low  programming  voltage.  The  analytical 
approximation  shows  relative  large  deviation  from  the  meas¬ 
ured  data  at  the  high-resistance  state. 

We  proposed  a  compact  model  to  simulate  the  transition 
region  motion  in  the  TiCT-TiCT  _x  memristor  based  on  the 
classic  ion  transportation  theory.  Our  model  is  validated  with 
the  measured  data  from  a  real  TiCL/TiCT  _x  memristor  device 
and  proved  capable  of  simulating  the  static  and  dynamic 


FIG.  3.  Model  validation  with  memristor  dynamic  switching  for  one  single 
cycle,  including  numerical  simulation  and  analytical  approximation. 


switching  properties  of  the  device  with  good  accuracy.  Our 
future  work  will  focus  on  modeling  the  impact  of  tempera¬ 
ture  and  3-D  model.  We  also  plan  to  publish  the  correspond¬ 
ing  Verilog-A  model  in  the  public  domain  in  the  near  future. 
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